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1 


This work addresses the question of the stability of stratified, spatially periodic shear 
flows at low Peclet number but high Reynolds number. This little-studied limit is 
motivated by astrophysical systems, where the Prandtl number is often very small. 
Furthermore, it can be studied using a reduced set of “low-Peclet-number equations” 
proposed by Lignieres [Astronomy & Astrophysics, 348, 933-939, 1999]. Through a 
linear stability analysis, we first determine the conditions for instability to infinites¬ 
imal perturbations. We formally extend Squire’s theorem to the low-Peclet-number 
equations, which shows that the first unstable mode is always two-dimensional. We 
then perform an energy stability analysis of the low-Peclet-number equations and 
prove that for a given value of the Reynolds number, above a critical strength of the 
stratification, any smooth periodic shear flow is stable to perturbations of arbitrary 
amplitude. In that parameter regime, the flow can only be laminar and turbulent 
mixing does not take place. Finding that the conditions for linear and energy stabil¬ 
ity are different, we thus identify a region in parameter space where finite-amplitude 
instabilities could exist. Using direct numerical simulations, we indeed find that 
the system is subject to such finite-amplitude instabilities. We determine numeri¬ 
cally how far into the linearly stable region of parameter space turbulence can be 
sustained. 
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I. INTRODUCTION 


The study of the onset of turbulence in stratified shear flows has a long history that dates 
back to Richardson 1 . He argued that the kinetic energy of turbulent eddies in a stratified 
shear flow can only decrease if N 2 /S 2 > 1, where N is the local buoyancy frequency, and 
S = \du/dz\ is the local shearing rate of the flow field u in the vertical direction e z . 
This criterion, derived simply from energetic arguments, is now commonly referred to as 
Richardson’s criterion, and the local ratio 


J(z) = 


N 2 (z) 

~S^z) 


( 1 ) 


is called the gradient Richardson number. The first linear stability analysis of a stratified 
shear flow is due to Taylor^, who considered both continuously and discretely varying strat¬ 
ification and shear profiles. This work, together with Goldsteinthen led to the derivation 
of the Taylor-Goldstein eigenvalue equation for the complex growth rate of two-dimensional 
infinitesimal disturbances in stratified shear flows. The solution of this equation for a given 
shear profile S(z) and stratification profile N(z) can be obtained either analytically in a few 
particular cases, or numerically in general. It wasn’t until much later, however, that the first 
general result on the stability of stratified shear flows was derived by Miles 4 and Howard 5 : 
a system is stable to infinitesimal perturbations provided J(z) is everywhere larger than 
1/4. As discussed by Howard and Maslowe^, this theorem should not be viewed as a refine¬ 
ment of Richardson’s argument (i.e. replacing 1 by 1/4), since the latter was specifically 
interested in determining when turbulence could be sustained, rather than triggered. In 
this sense, Richardson’s original argument should be viewed more as a nonlinear stability 
criterion than a linear one. 

These results were obtained in the limit of vanishing viscosity and diffusivity. For ther¬ 
mally stratified flows, however, thermal diffusion can have a significant influence on the 
development of shear instabilities by damping the buoyancy restoring force. This effect was 
first studied by Townsend 1,1 in the context of atmospheric flows. He showed that the thermal 
adjustment of the fluid parcel to its surroundings, by radiative heating and cooling or by 
thermal conduction, always acts to destabilize the flow and increases the critical Richard¬ 
son number for linear stability, J C rit, by a factor inversely proportional to the product of 
the shearing rate S with the cooling time t coo \ (this product is a local Peclet number for 
the flow), so J crit ~ (tcooiS') -1 . Viscosity, meanwhile, has a generally stabilizing influence. 



Zahn® emphasized the importance of these results for stellar astrophysics: in stellar interiors 
where the Prandtl number is typically very small (Pr ~ 10~ 8 — 10 -5 ) high Reynolds number 
flows can also have a low Peclet number, or in other words, thermally diffusive shear flows 
exist when viscosity is nevertheless small enough not to suppress the development of the 
instability. This combination is ideal for shear instabilities, and is specific to astrophysical 
systems - it cannot happen for most geophysical flows where the Prandtl number is usually 
of order unity or larger. 

Applying Townsend 7 ^s results to shear-induced turbulence in stellar interiors, Zahn 9 fur¬ 
ther argued that the relevant cooling timescale is the radiative timescale based on the size 
l of turbulent eddies, namely f coo i = l 2 / Kt where Kt is the thermal cliffusivity. He then 
proposed to take for l the smallest length scale for which viscosity is still negligible, that is, 
one for which the turbulent Reynolds number Re; = Sl 2 fv = St coo i/ Pr ~ Re crit (where v 
is the kinematic viscosity), where Re cr i t is a constant that he estimates to be around 10 3 . 
This would imply J crit ~ (SAooi)” 1 ~ Re“j t P r_1 , or hi other words, J cr i t Pr ~ Re“( t ~ 10~ 3 . 
Zahn®s argument, as in the case of Richardson 1 ^ original argument, should be viewed as 
a nonlinear stability criterion rather than a linear one, since it relics on the presence of 
pre-existing turbulent eddies. 

Zahn’s work had an enormous impact in the field of stellar evolution. While the standard 
Richardson criterion is far too stringent to allow for the development of shear instabilities 
in the absence of thermal diffusion (the typical Richardson number being much larger than 
one even in the strongest known stellar shear layers), its relaxation allows for the possibility 
of much-needed mixing in stellar evolution theory. Indeed, models without any form of 
turbulent mixing in stably stratified regions are not able to account for observations. As 
reviewed by Pinsonneault 10 the problem is particularly acute when it comes to explaining 
the surface chemical abundances of light elements such as lithium and beryllium, as well 
as products and by-products of nuclear reactions such as helium, carbon, nitrogen and 
oxygen. More recently, further indication of the need for turbulent mixing was revealed 
by asteroseisomolgy thanks to the Kepler mission. Measurements of the internal rotation 
rate of red giant star^® 1 are inconsistent with evolution models in which turbulent angular- 
momentum transport in stably stratified regions is neglected. In both cases, therefore, 
efficient chemical transport and angular-momentum transport by shear instabilities could 
be the key to resolving these problems - the question remains, however, of whether these 


instabilities are indeed triggered, and how efficient mixing is. 


In the limit of low-Peclet numbers (i.e., high thermal diffusivity), the temperature fluctu¬ 
ations are slaved to the vertical velocity. The corresponding quasi-static approximation was 
originally introduced to study low-Prandtl-number thermal convection 1213 . In the context 
of stably-stratified systems, the quasi-static approximation was introduced only recently by 
Lignieres 14 (see Section IIC for more detail). He showed that the standard Boussinesq equa¬ 
tions can be replaced by a reduced model that is valid in the asymptotic low-Peclet-number 
limit, and that this model only depends on two parameters: the Reynolds number, and the 
product of the Richardson number with the Peclet number. As a result, the linear stability 
properties of the system depend on the product PegJ (where Peg = SL 2 s /kt with Lg being 
a characteristic vertical length scale of the laminar flow) rather than on each parameter 
individually. Since Peg is small by assumption, shear-induced turbulence can be expected 
even if the Richardson number is much larger than one. 


By contrast with linear theory, very little is known to date about the stability of stratified 
shear flows to finite-amplitude perturbations when viscosity and thermal diffusivity are both 
taken into account. It is yet a question of crucial importance in stellar astrophysics, since the 
presence or absence of vertical mixing can strongly affect model predictions. We therefore set 
out in this work to characterize the domain of instability of low-Peclet-number shear flows 
to finite-amplitude perturbations. For simplicity, we consider a specific shear profile that is 
periodic in the vertical direction. We present the model in Section [TT[ and briefly discuss 
the low-Peclet-number asymptotic equations proposed by Lignieres^ 14 . In Sections III and 


IV, we study the linear and nonlinear stability of the system respectively, and contrast the 


results in the low-Peclet-number approximation to those obtained starting from the full set 
of primitive equations. As we shall demonstrate using an energy stability analysis of the low- 
Peclet number equations, smooth periodic shear flows are stable to perturbations of arbitrary 
amplitude for sufficiently large Richardson number, for a given value of the Reynolds number. 
In Section [VJ we turn to direct numerical simulations to study the transition to turbulence 
via linear instabilities and finite-amplitude instabilities. We summarize our results and 
conclude in Section ED 







II. 


THE MODEL 


A. Model setup 



FIG. 1. Model setup: a horizontal shear flow is driven by a body-force. The background stratifi¬ 
cation is linear, and the temperature and velocity fluctuations are periodic in the three directions. 

Since our intention is to study the energy stability properties of stratified shear flows, it 
is crucial to start with a model where the mechanism driving the shear is explicit, which 
guarantees a well-defined energy budget. Two options are available: boundary-forcing, and 
body-forcing. Having potential applications to stellar astrophysics in mind, we prefer the 
latter in order to avoid boundary layer dynamics near solid walls, which are rarely present 
in stars. 

A simple and numerically-efficient way of studying body-forced, stratified shear flows is 
to consider a Boussinesq systemf 1 ^, where the forcing and all the perturbations are triply- 
periodic, and where the background density is linearly stratified-® (see Figure [T|. The model 
equations describing such a system are: 

du 1 1 

— —b u • Vm = - Vp + agTe z + v\/ 2 u -\ - F , (2) 

dt p 0 p 0 

V • u = 0 , (3) 

dT 

— +u-VT + wT 0z = k t V 2 T , (4) 

where u = ( u , v, w) is the triply-periodic velocity field, p and T are the triply-periodic 
pressure and temperature perturbations, po is the mean density of the region considered, a 
is the coefficient of thermal expansion, g is gravity, v and kt are the viscosity and thermal 
diffusivity (respectively). The quantities po, a , g, z/, k are all assumed to be constant, 
as in the standard Boussinesq approximation. The use of the latter is justified as long 












as the vertical height of the domain is much smaller than a density scalehcight. Finally, 
we assume that there is a constant background temperature gradient. 1 ' T 02 , and that all 
thermodynamical and dynamical perturbations have zero mean in the domain. 

The applied force F should be triply-periodic as well. A natural candidate is a sinusoidal 
forcing, thus driving a Kolmogorov flow in the laminar regime. In what follows, we assume 
that F is of the form 

F = F 0 sin(kz)e x , (5) 


which defines a typical lcngthsca.le k h In the steady laminar regime, this force generates a 
sinusoidal shear flow along the ^-direction, 


u L 


Fq 

p 0 uk 2 


sin (kz) e x . 


( 6 ) 


Note that while the present paper deals mostly with Kolmogorov forcing, the energy stability 
of arbitrary smooth velocity profiles is discussed in section |IVB 2 


B. Non-dimensionalization and model parameters 


We non-dimensionalize the equations using the amplitude of the laminar solution, ^° fc2 , 
as a velocity scale. We also use the spatial scale of the laminar solution, k~ 1 , as the unit 
length scale. This then defines the timescale With this choice of units, the equations 
©-0 become 


du 

~dt 

V • u — 0 , 
dT 

m +uVT 


——b u • Vm = — Vp + RiTe 2 


w = i^ v2r ■ 


yJ-V 2 (u - sin(^)e x ) , 
Ke 


( 7 ) 

( 8 ) 

(9) 


where 


agT 0z plv 2 k 2 


O _ 10 TT — ^ -n _ 1 ^ /im 

° “ p 0 u 2 k 3 ’ 1 “ F 2 ’ e “ pouk 3 Kt • ( ) 

The laminar solution (|6]) now takes the dimensionless form ui = sin(z) e x . Provided the 
system remains in the vicinity of this laminar solution, Re, Pe and Ri are the usual Reynolds, 
Peclet and Richardson numbers based on the typical velocity of the flow. 










C. Low Peclet number approximation 


When a field diffuses on a timescale much shorter than the advective time, it enters a 
quasi-static regime where the source term and diffusive term instantaneously balance. Such 
a quasi-static regime has been used for decades in the context of magneto-hydrodynamics 
of liquid metals: at low magnetic Reynolds number, the induced magnetic field is slaved to 
the velocity field 18 . The equivalent approximation for flows of low Prandtl number fluids 
was originally introduced in the context of thermal convection 1213 . Rather surprisingly, 
this quasi-static approximation appeared only much more recently in the context of stably- 
stratified flows. 

Lignieres 14 proposed that in the low-Peclet-number limit the governing equations ([7])-([9]) 
can be approximated by the reduced set of so-called “low-Peclet-number” equations (LPN 
equations hereafter) 


du 

——b u • Vm = 
at 

V • u = 0 , 

w = v^ T - 


-Vp + RiTe^ + —V 2 u + F , 
Re 


( 11 ) 

( 12 ) 

(13) 


to zeroth and first order in Pe. To derive equation (13), one can assume a regular expansion 
of T in Pe, as in T = To + PeT) + 0(Pe 2 ), and further assume that the velocity field is 
of order unity. At the lowest order, the temperature equation yields V 2 Xo = 0 which then 
implies T 0 = 0 given the applied boundary conditions. The next order then yields the 
quasi-static balance w = V 2 7\ ~ Pe _1 V 2 T as required. It is worth mentioning that in 
Lignieres®s original work the velocity is scaled with its dimensional r.m.s. value u rms , so 
it is Pe rms = u rms /kKT, rather than Pe, that has to be small for the LPN equations to be 
valid. In Sections H and |IV[ we shall study the linear and energy stability of a laminar 
flow for which the r.m.s. velocity is of the same order as the flow amplitude. In that case 
Pe ~ Pe rms and we expect the LPN equations to be valid whenever Pe is small. In Section 
[VJ however, we shall see that numerical simulations of turbulent shear flows that have large 
Pe can still be well-described by the LPN equations as long as Pe rms is small. 

It is also worth noting that the LPN equations only allow for temperature fluctuations 
T that have a zero horizontal mean, by contrast with the full equations. As a result, the 
horizontal mean of the full temperature field (background plus perturbations) is necessarily 





linear in z. To see this, we take the horizontal average of the thermal equation, which, 
assuming that there is no vertical mean flow (which can be guaranteed by making sure the 
initial conditions do not have one), results in 


d 2 (T) h 
dz 2 


0 , 


(14) 


in the LPN equations, where (T)h is the horizontal average of T. The only solution of 
this equation which satisfies periodicity is the constant solution; further requiring that the 
volume-average of T be zero then implies that (T)h = 0. By contrast, taking the horizontal 
average of the standard temperature equation under the same assumptions would result in 


d(T) h 

dt 


+ 


d 

dz 


{wT) h = 


1 d 2 (T) h 
Pe dz 2 


(15) 


which has solutions with non-zero (T)^. These solutions can, for instance, develop into 
density staircases under the right circumstances 19 20 . The latter are however prohibited 
in the LPN equations. Among other effects, this rules out the development of Holmboe 
modern, and may explain why it is possible to get simple energy stability results in the 
LPN limit but not for the full equations. 

Finally, combining the momentum and the thermal energy equations yields 
du _ _ o 1 ,, 

— — hit' Vm = — Vp + RiPeV ~we z + —- V 2 u + F , (16) 


which formally shows that the Richardson number is no longer the relevant parameter of the 
system, but that RiPe is. The work of Lignieres® thus puts the arguments of Townsend 7 
and Zahn 9 discussed in Section [I] on a firm theoretical footing. Lignieres 14 and Lignieres, 
Califano, and Mangeney 22 verified that the LPN equations correctly account for the linear 
stability properties of various systems in the low-Peclet-number limit. Prat and Lignieres® 
later also verified that they correctly reproduced the low-Peclet dynamics of their 3D non¬ 
linear simulations. In this paper, we continue to verify the validity of the LPN equations 
through stability analyses and nonlinear simulations. 







III. LINEAR STABILITY OF A PERIODIC KOLMOGOROV FLOW 


We first focus on the stability of the laminar solution to infinitesimal perturbations. We 
solve the linearized versions of equations 0-0, 

f)u' 1 

— + u L -Vu' + u'- Vu L = - Vp + RiT'e 2 + —V V , (17) 

at Re 

V • v! = 0 , (18) 

% + u L ■ vr + w‘ = ^ w , (19) 

where u' and T’ are infinitesimal perturbations to the linearly stratified background shear 
flow ul = singes;. 


A. Squire’s transformation 


The linear stability of the unstratified Kolmogorov flow was first investigated in detail 
by Beaumont. 24 . Squire’s theorem 25 states that the first unstable mode as the Reynolds 
number increases is a (^-independent) 2D mode. This strong result implies that one can 
focus on 2D perturbations to determine the stability threshold of the system. Such a 2D 
analysis is much simpler and computationally less expensive than a 3D one. Beaumont 21 
found that 2D flows are unstable only for Re > a/2- 

The linear stability of the stratified Kolmogorov flow to 2D perturbations was studied 
in detail by Balmforth and Young 26 . The 2D case can be made more generally relevant 
by noting that Squire’s transformation 25 for the viscous unstratified case can be extended 
to the stratified case with thermal diffusion to argue that the linear stability of any 3D 
mode can equivalently be studied by considering that of a 2D mode at lower or equal 
Reynolds and Peclet numbers, and higher or equal Richardson number. This result, which 
was summarily discussed by Yih 27 and clarified by Smyth, Klaassen, and Peltier 28 and 
Smyth and Peltier^®, states that the growth rate A 3 of the 3D normal mode qs(x, y , z, t ) = 
qs(z)exp(ilx + imy + A 3 t) at parameters (Re, Pe, Ri) is related to that of the 2D normal 
mode q 2 (x, y , z, t ) = q-iiz) exp (iLx + A 2 1) at suitably rescaled parameters via 

A 3 =/(/,m;Re,Pe,Ri) = ^A 2 =-^/^L,0;-^-Re,-^-Pe, ^-rA , (20) 


where L = y/l 2 + m 2 . One can apply the same method to the LPN equations, and the 


result can readily be deduced from (20). Indeed, the transformation (20) is valid for any 







Peclet number, so it remains valid for low Peclet numbers. In this limit, we saw that only 
the product RiPe is a relevant parameter: as a consequence, one can replace the last two 


arguments of the function / in (20) by the product of the two. The low-Peclet version of 


Squire’s transformation therefore gives 


A 3 = film: Re, RiPe) = 2 A 2 = )-/ (i, 0; )-Re, yRiPe 


( 21 ) 


This relationship between the growth rates of 2D and 3D modes has important impli¬ 
cations for the marginal linear stability surface. In order to find the latter in 2D, we first 
maximize the real part of f(L , 0; Re, Pe, Ri) over all possible values of L , yielding the function 
5^ (Re, Pe, Ri) which returns the growth rate of the fastest growing mode for each parameter 
set (Re, Pe,Ri). The marginal linear stability surface is then defined by S 2 = 0. Similarly, 
the marginal linear stability surface for 3D perturbations is obtained by constructing the 
function S^Re, Pe, Ri) = max i]m Re[f(l, m; Re, Pe, Ri)] and setting S 3 = 0. If the functions 
S 2 and S 3 are the same, then so are the surfaces S 2 = 0 and S 3 = 0 , which implies in turn 
that the first modes to be destabilized are the 2D modes. This is the case for instance in 
the limit where stratification is negligible (see above). 

In general, the only way to determine whether S 2 = S 3 for a given linear stability problem 
is to construct these functions by brute force, using their original definition as the growth 
rates of the fastest growing modes. While this is not too time-consuming in 2D, it can 
become computationally expensive in 3D. However in this particular problem, since the 
growth rates of 2D and 3D modes are related, we also have: 


S 3 (Re, Pe, Ri) = max yS' 2 (yRe, yPe, Ri/y 2 ) , ( 22 ) 

xe[o.i] 

where y = \l/L\. A similar relationship applies for the LPN equations: 


S^(Re, RiPe) = max yS'AyRe, RiPe/y) . 
xe[o,i] 


(23) 


Note that it is easier to construct S 3 from equations ( 22 ) or (23) than to do so directly. 

Whether S 2 = S 3 or not then simply depends on the properties of S 2 . It is quite easy to 
find sufficient conditions that guarantee S 2 = S 3 . For instance, in the case of the standard 
equations, if S 2 is a strictly increasing function of both Re and Pe, and a strictly decreasing 
function of Ri, then the maximum over all possible values of y is achieved for y = 1, which 
ensures that S 2 = S 3 . For the LPN equations, it is sufficient to show that S 2 is a strictly 





increasing function of Re and a strictly decreasing function of RiPe. In what follows, we 
therefore first study the stability of 2D modes, and then use these results to conclude on the 
stability of the system to 3D modes. 


B. Linear stability analysis using Floquet theory 


We use a stream function to describe divergence-free 2D perturbations, 


u = u L + V x pip'e y ) , (24) 

where 0' is the infinitesimal perturbation. The linearized equations (|I7|)-(|T9| become 

|(W)+ *.(*) (|(VV+^)) = + ivy , 


dT . , dT' <90' 1 

-+ sin z — + w- = —W 


(25) 


dt dx dx Pe 

This set of PDEs for T' and ip' has coefficients that are independent of t and x, but periodic 
in z. Normal modes for this system are of the form 


d(x,z,t)=e iL *+ xt q(z) , 


(26) 


where q' is either T' or 0', and L is real. Using Floquet theory, we then seek solutions for q 
given by 

N 

q(z) = e‘“ V , (27) 

n=—N 

where a is real, to satisfy the general periodicity of the system. Substituting this ansatz into 
the previous equations, we obtain an algebraic system for the and T n : 

—A((a + n ) 2 + L 2 )ip n + — [(1 — (a + n — l) 2 — L 2 )ip n -i — (1 — (a + n + l) 2 — L 2 )0 n+ i] 

= iK\LT n + z— ((a + n) 2 + L 2 ) 2 0 n , 

Re 

ATn + [T n ~ i — T n+ 1] + iLip n = n y L 2 )T n , (28) 

for n = —N...N. This can be cast as the linear eigenproblem, 


M(L; a; Re, Pe, Ri)£c = \x , (29) 

where x = {pJ-N, t/hv, T-jv, •••, Rv}, which can be solved for the complex growth rate 
A. The real part of the latter can then be maximized over all possible values of a and L 





for given system parameters (Re, Pe, Ri) to determine the temporal behavior and spatial 
structure of the most rapidly growing mode of the shear instability, or in other words, to 
construct S 2 (see previous section). In both unstratified and stratified cases studied so far, 
the first unstable modes at the instability threshold have the same periodicity in z as that of 
the background shear, so that a = G 212631 . We verified that this is indeed the case here as 
well. In what follows, we therefore restrict the presentation of our results to the case a = 0. 

The equivalent problem for the LPN equations is given by 

— A((a + n)“ + L 2 )Yn + 7^ [(1 — (a + u — l) 2 — L 2 )Yn-i — (1 — (a + u + l) 2 — L 2 )Yn+i] 

= + Y + + k ((a + nf + i2)V " ’ (30) 

which can be cast as 

M'(L; a; Re, RiPe)?/ = A y , (31) 

where y = {Y-n, ■ This time, the fastest growing mode only depends on two system 

parameters, namely Re and the product RiPe. Again, we restrict the following analysis to 
the case a = 0. 

Various aspects of the marginal stability surface S^Re, Pe, Ri) = 0 for 2D modes are pre¬ 
sented in Figure [2] Figure^ shows the critical Reynolds number as a function of Ri for the 
standard equations, for various values of the Prandtl number (Pr = Pe/Re). The evolution 
of the shape of these curves as Pr increases is not a priori easy to identify nor explain. How¬ 
ever, an obvious result is the existence of unstable modes for reasonably large values of the 
Richardson number when the Prandtl number is low. This can easily be understood in the 
light of the work of Townsend 7 (see also Gage and Millei®, Jones 33 ,Lignieres 14 ,Lignieres, 
Califano, and Mangeney 22 for instance), who showed that stratified shear instabilities can 
exist beyond the standard Richardson criterion when thermal diffusion is important. Since 
thermal diffusion increases as Pe decreases, and since Pe = Pr Re, one can naturally expect 
unstable modes at high Richardson number for fixed Re and low enough Pr (or vice-versa). 

Following Lignieres, Califano, and Mangeney 2 ^, we now show in Figure^ the same data 
plotted against RiPe, and add the marginal stability curve for the LPN equations. The 
interpretation of the results is now much clearer. We first see that the LPN equations are 
indeed a good approximation to the full equations when the Peclet number is small (low 
PrRe). I 11 the unstratified limit Ri —> 0, we find that marginal stability is indeed achieved 







FIG. 2. Left: Marginal stability curves for 2D modes in the form of Re vs. Ri for various values of 
the Prandtl number: Pr = 10 -3 (green dashed line), Pr = 10 -2 (blue small dashed line), Pr = 10 -1 
(purple dotted line), Pr = 1 (cyan long dash - dotted line), and Pr = 10 (brown short dash - dotted 
line). The system is unstable in the area above and to the left of the curves. Right: the same 
data plotted against RiPe instead. The red solid line is the marginal stability curve for the LPN 
equations. The Pr = 1CU 3 and Pr = 1CP 2 curves nearly overlap with it for Re < 100, which is 
consistent with the fact that Pe < 1 for these values of the Prandtl number. In all cases, we have 
truncated the Fourier expansion of i\)' and T' to N = 20 to create these curves. This choice of N 
was made after successful convergence tests. 

for Re = y/2, as expected 24 . We also see that, for the low-Peclet equations, the threshold for 
linear stability (RiPe)L above which the flow is linearly stable becomes independent of the 
Reynolds number for large enough Re. The asymptotic value can be estimated numerically, 
and is roughly (RiPe^Rg-^oo — 0.25. The fact that the critical RiPe for linear stability is 
independent of the Reynolds number for large enough Re shows that the inviscid limit is 
a regular limit of this problem. This is, however, in contrast with the findings of Jones^ 3 
and Lignieres, Califano, and Mangeney 22 for the tanh shear layer. In both cases, they find 
that (RiPe)i oc Re for large Re (albeit using a fairly limited survey of parameter space). 
Using the LPN equations, Lignieres, Califano, and Mangeney 22 also found that there is no 








stability threshold in the inviscid limit, that is, unstable modes exist for all values of RiPe. 
The reason for the stark difference between our results and theirs remains to be determined, 
but could be attributed either to the nature of the boundary conditions used (periodic vs. 
non-periodic), or to the fact that a sinusoidal velocity profile has shear of both signs while 
a tanh velocity profile only has shear of one sign. 

We now discuss linear stability to 3D perturbations using Squire’s theorem. For the LPN 
equations, we find that the function £2 (Re, RiPe) is indeed a strictly increasing function of 
Re, and a strictly decreasing function of RiPe, which implies that the marginal stability of 
2D modes is also that of 3D modes (see the previous section). For larger Prandtl number, 
however, we can immediately see from its null contour that £2 (Re, Pe, Ri) is no longer a 
monotonic function of RiPe which strongly suggests that 3D modes could be the first ones 
to destabilize the system. Whether this is indeed the case is beyond the scope of this paper, 
since it belongs to the high-Peclet-number regime. However, this result would be consistent 
with the work of Smyth and Peltier 29 [ who found that 3D modes can be the first ones to be 
unstable for parallel stratified shear flows which have a tanh profile, albeit in some relatively 
small region of parameter space. 

In preparation for our 3D simulations (see Section |VJ), we are also interested in the spatial 
structure of the first modes to be destabilized, since the computational domain size must be 
chosen to be large enough to contain them in order to avoid spurious results. Based on the 
previous results, we now limit our study entirely to the 2D modes. The range of unstable 2D 
modes for which marginal stability is achieved is shown in Figure [3j for both the standard 
equations and for the LPN equations. Again, we see that the results obtained using the 
LPN equations correctly approximate those obtained using the standard equations at low 
Peclet number. In all cases, we find that the first mode to be destabilized has L e [0,1], 
i.e. its horizontal wavelength is larger than the shear lengthscale. For this reason, in the 
numerical simulations of Section [V] we use a reasonably long domain size, that can fit at 
least two wavelengths of the most unstable mode. 

Finally, Figure [3] also sheds light on the actual source of the non-monotonicity of the 2D 
linear stability curves seen in Figure[2]at Pr = 1 and above. Indeed, it reveals a new unstable 
region for low horizontal wavenumber modes, which appears here for Pe = 10 (Pr = 0.1 in 
this figure). These modes have a growth rate A with non-zero imaginary part, by contrast 
with the standard shearing modes whose growth rates are reaP 4 


100 



FIG. 3. Marginal stability curves for 2D modes, for Re = 10 2 , in the form of RiPe vs. L, where L 
is the horizontal wavenumber of the first unstable mode. The system is unstable in the area within 
the curves, which therefore corresponds to the range of unstable modes for a given Richardson- 
Peclet number. The solid red line was obtained using the LPN equations. The green dashed line 
is for Pe = 0.1 (Pr = 10 -3 ), the blue short-dashed line for Pe = 1 (Pr = 10 2 ), the purple dotted 
line for Pe = 10 (Pr = 10” 1 ) and the cyan dot-dashed line for Pe = 100 (Pr = 1). As in Figure [ 2 J 
we have truncated the Fourier expansion of if and T' to N = 20 to create these curves. 

IV. ENERGY STABILITY 

A. Energy stability for stratified shear flows: general ideas 

Linear stability only provides information on the stability of a shear flow to infinitesimal 
perturbations. Energy stability is a much stronger form of stability 30 35 36 : when a system 
is energy stable (also called absolutely stable), perturbations of arbitrarily large amplitudes 
decay at least exponentially in time, and the laminar flow is the only attractor of the system. 
Energy stability is thus a sufficient condition for the system to be stable to perturbations 
of arbitrary amplitude, whereas linear instability is a sufficient condition for the system to 
be unstable. Often the linear stability limit and the energy stability limit do not coincide 
in parameter space, an extreme example being the unstratified plane Couette flow, which 
has a finite threshold Reynolds number for energy stability, but is linearly stable up to 







infinite Reynolds number: in the region between the two, the system is stable to infinitesimal 
perturbations, but it may be unstable to perturbations of large amplitude, i.e., it may exhibit 
finite-amplitude instabilities. 

With the goal of further studying the stabilizing effect of background stratification, we 
now derive an energy stability criterion for forced stratified shear flows. We ask the following 
question: for a given amplitude of the force, is there a critical strength of the stratification 
above which the laminar solution is the only attractor of the system? 

Again, we insert the decomposition u(x,t) = u L (z ) + u'(x,t) in (J7])-(J9]). However, we do 
not assume that u' is small, u' and T' then satisfy 


——b ul • Vw' + u • Vul + u • Vu 7 

at 

V • u' — 0 , 


dV 

~dt 


+ (u L + u ') ■ VT 7 + w' 
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V 2 T 7 
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VV 


(32) 

(33) 

(34) 


in the case of the standard equations, and 

3u f 1 

——b u L • Vm' + u • + u • Vm' = — V» + RiPeV _ 2 w , e 2 + —V 2 u' , (35) 

ot Re 

V • v! = 0 , (36) 


for the LPN equations. 

An energy equation for the perturbations can be obtained by dotting the momentum 
equation with u', adding it to yT 7 times the temperature equation (where the only constraint 
on 7 is that it should be a positive scalar), and integrating the result over the domain under 
consideration. Using the periodicity of the solutions, together with the incompressibility 
condition greatly simplifies the resulting expressions, which reduce to 

d_ 
dt 

= (37) 


)( U ' 2 ) + ! (T“> 


= (Ri - 7 ){w'T') - (S L w'v!) - jhflVuf) - j^dVTf) , 


for the standard equations, where (•) denotes a volume integral, and Sl(z ) = ^Ml(^) = 
cos(z) denotes the local vertical shear of the laminar solution. Similarly, for the LPN 
equations we get 


d_ 

dt 



RiPe(ta 7 V~ 2 ta 7 ) - (S L w'u') 
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Re 
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This defines the two functionals TL- [u',T'] and T-Llpn [w 7 ], which are both quadratic forms. 
The task at hand is to determine the region of parameter space {Re, Pe, Ri} or {Re, RiPe} 
where these quadratic forms are negative definite, i.e., where they are strictly negative for 
any possible input fields u' (and T'). In this region of parameter space the system is energy 


stable: the right-hand-side of (37) or (38) is strictly negative and the perturbation decays in 
time, regardless of its initial amplitude. On the basis of mass and momentum conservation, 
the only constraints that we place on u' is that it is divergence-free and has a vanishing 
average over the whole domain. 


Comparing 'H 7 [u',T'] and / Hlpn[ u '], we readily see that the task of proving energy 
stability for stratified shear flows is more involved in the case of the full set of equations 
than in the case of the LPN equations. We therefore focus on the latter, for which we are 
able to obtain interesting and generic results on the stability of stratified shear flows. 


B. Energy Stability in the low Peclet number limit 

Focussing on the LPN equations, we first compute bounds on the location of the energy 
stability curve in the (RiPe, Re) plane. These bounds prove useful, because they correctly 
describe the scaling behavior of the energy stability limit for large Reynolds number. As we 
shall see they also validate our numerical results, and provide a simple analytical approxi¬ 
mation to the high Reynolds number limit. 


1. Lower bound 

While the true energy stability boundary can only be obtained by ensuring that H LPN < 0 
for all possible perturbations u' and T', one may also ask the question of when a subset of 
perturbations is energy stable or unstable. If the subset is unstable, then we know that the 
system overall is not energy stable either. The critical RiPe for energy stability for that 
subset is therefore a lower bound (called (RiPe)< hereafter) on the true energy stability 
boundary. 




We now restrict our attention to perturbations of the following form: 


u' = B cos (JCy ), 

v' = — sin(/C?/) sin ( 2 ), 

1C 

vJ = cos(z) cos(ICy). 


(39) 

(40) 

(41) 


where /C is the wave number of the perturbation in the y direction, and B is a free parameter. 
One can check that such perturbations are divergence-free. 


We insert (39)-(41) into the quadratic form (38), recalling that Sl = cos(z), to obtain 

B 1 r IC 2 B 2 
J-Re 


Blpn 
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L x LyL z 
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1 + K 1 


(42) 


This expression is a quadratic polynomial in B. As long as its discriminant is negative, 


the polynomial is negative and perturbations of the form (39)-(41) decay exponentially. 
However, when the discriminant is positive, there will be values of B corresponding to 


growing perturbations. The threshold for energy stability of perturbations of the form (39)- 


(43) 


(41) is therefore attained when the discriminant vanishes, which gives 

c 2 +1„ (K 2 + l) 2 (l+ i) 

(RlPe) < = ^ Re -' 

For large enough Reynolds number, this value is maximum for the smallest value of /C that 
is compatible with the boundary conditions, namely /C = 2n/L y . The corresponding lower 
bound on the energy stability limit of the system is 


(RiPe)< = 


R eL 2 v 
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For a given system size, the high-Re asymptotic behavior of the lower-bound is 


2 3 L l 

(RiPe)< ~ y 


32t r 2 


4tt 2 

-jY + 1 ) Re , 


(45) 


which shows that the Richardson-Peclet number needs to be at least of the order of Re to 
have energy stability. This bound also indicates a strong dependence of the energy stability 
limit on the size of the domain. Indeed, as the transverse size L y of the domain increases, 


expression (45) grows as L 2 : larger domains allow for perturbations that are very weakly 
damped by viscosity. 

The lower bound is plotted in Figure [4] for a domain of size 107r x 2n x 27T, for which 

Re 8 
= T ~ RW 


(RiPe)< 
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2. Upper bound 


Upper bounds on the energy stability limit can be obtained using rigorous estimates of 
the three terms in T-Llpn■ In this subsection, we do not restrict attention to shear flows of 
the Kolmogorov type. Instead, we consider any smooth shear flow ul(z) along x that is 
27r-periodic in z. We still use the height of the domain as the characteristic length scale, 
and consider a domain of size L x x L y x 27r with periodic boundary conditions. We assume 
that some forcing function with amplitude F 0 sustains the laminar flow. The dimensionless 
profile has an amplitude of order unity. We prove that any such laminar flow is energy stable 
provided the Richardson-Peclet number is large enough. 

To simplify notations and avoid dealing with the inverse Laplacian operator, we introduce 
9 = T'/Pe, such that w' = V 2 R With these notations, the quadratic functional reads 

U LPN [«'] = —RiPe(| Vd| 2 ) - (^VV 2 0> - ^(|Vu'| 2 >. (47) 


Let T be the second term of this functional. After integration by parts, 
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Using classical inequalities, we bound this term according to 
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iwn, ( 49 ) 


where we have used respectively Holder’s inequality to get the first line, Cauchy-Schwartz 
inequality to get the second one, and Young’s inequality to get the final expression (see 
Doering, Spiegel, and Worthing® 7 for another example of the use of these inequalities in a 
fluid dynamics context). We now wish to express (| , u / | 2 ) in terms of (| ’Vu'\ 2 ). To wit, we use 
Poincare’s inequality 3 ': the divergence-free constraint implies that u’ has vanishing Fourier 









































amplitude on modes with vanishing wave numbers in both the y and z directions, hence 


<I«T) < ^max{Lj;4ir 2 } (|Vt/| 2 ), 


(50) 


where the arguments of the max are the maximum allowed values for the squared wavelengths 


in the y and z directions. Inserting this inequality into (49), together with (|W| 2 ) < 
(IV-u'I 2 ), leads to 
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As discussed at the beginning of this subsection, the non-dimensionalization is such that 

(52) 

(53) 


sup 


d 2 u 2 


= Cl , 
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= c 2 , 


where Ci and C 2 are constants of order unity that depend on the shape of the laminar profile 


only (and specifically not on Re, RiPe, L yi etc). Combining these expressions with (51) 
leads to 


7-Llpn [w'j < <|Vu'| ) 


2 RiPe 
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Re 


(54) 


hence 1-Llpn is a negative quadratic form if the expression inside the square brackets is 
negative, i.e., if 


RiPe > (RiPe)> = ( ci + c 2 max <( 1 \ ) ^ . 


(55) 


Because we have used rough but rigorous estimates of the different terms of the quadratic 
functional, (RiPe)> is an upper bound on the actual energy stability limit of the system. It 
proves that any smooth laminar velocity profile is absolutely stable provided the stratification 
is strong enough. This has important implications, showing in particular that for fixed 
Reynolds number and strong enough stratification such a shear flow does not induce any 
turbulent mixing! Note, however, that since this result is obtained for the LPN equations, 
it is formally only valid for perturbations that have a low Peclet number. This can be 
done mathematically by taking the asymptotic limit of the equations for Pe —)■ 0 before 
considering perturbations of arbitrary amplitude. In practice, however, our result does not 
rule out the possibility of instability for perturbations that have a high Peclet number. A 






















simple example would be perturbations that locally reduce or eliminate the horizontally- 
averaged vertical stratification the domain: such perturbations are not allowed in the LPN 
equations, but they are allowed in the full set of equations at very low Peclet number, where 
they might indeed allow for sustained turbulent solutions localized in the mixed layer. 


For domains with large extent in the ^-direction, the sufficient condition (55) for energy 


stability becomes approximately RiPe > L 2 Re. In the particular case of the Kolmogorov 


velocity profile, we can compare this upper bound to the lower bound (45): both bounds 


scale as L 2 Re for large Reynolds number, which indicates unambiguously that the actual 
critical RiPe for energy stability obeys the same scaling. This is illustrated in Figure [4| 
where we plot the upper bound for a Kolmogorov flow in a domain of size 107T x 27rx 27r. 
For such a Kolmogorov flow, C\ — c? — 1 and the bound becomes 


(RiPe)> = Re. 


(56) 


3. Energy stability boundary using Euler-Lagrange equations 

To determine the true energy stability limit of the LPN system, we now consider the 
variational problem associated with the extremization of the quadratic functional. This 
gives a set of Euler-Lagrange equations, which can be solved numerically to obtain the 
critical value of RiPe for absolute stability, called (RiPe)# hereafter. 

Starting from LLlpn [w 7 ], we separate the viscous dissipation term from the other two, as 

U LPN [u'\ = RiP'e(u/V~V} - {S L w'u') - V , (57) 

where T> = ^(iVit'l 2 ) is positive definite for non-trivial flows. We then ask the following 
question: for fixed viscous dissipation rate V = D 0 , for what values of RiPe and Re is 
LLlpn M < 0 for all incompressible velocity fields u'l As we shall see, the selected value 
of Do is irrelevant as it merely serves as a general normalization 38 of u!. In order to answer 
this question, it is sufficient to maximize RiPe(u/V _2 u/) — ( Slw'u') over all possible incom¬ 
pressible flows u ', and find out for what values of RiPe and Re this maximum is smaller 
than Do- The problem thus reduces to an Euler-Lagrange optimization problem. 

Using the notation 6 = T'/Pe as in the previous section, we construct the following 




Lagrangian: 


C [u'j = RiPe(u/#) — (Sz J w , u , ) + (7Ti(x, y , z) V- u') — 7r 2 ( 


I Vu' 


Re 


D 0 ) + {x 3 (x,y, z)(w'-V 2 6)) , 

(58) 

where the Lagrange multiplier function iri(x,y,z) enforces incompressibility at every point, 


iT 3 (x,y,z) enforces equation (13) at every point, and the constant multiplier 7r 2 enforces 


T> = D 0 globally 39 . Note that we go back here to using the field 9 merely to avoid dealing 
with inverse Laplacian operators in the variational problem; it is also possible to work 
through the following derivation without doing it. 

The maximizing perturbation field u' has to solve the Euler-Lagrange equations: three 
of them (arising from the derivatives of C with respect to 7Ti, 7r2 and 7 t 3 ) simply recover the 
constraints, and the other four (arising from the derivatives of C with respect to u', v', w' 
and 6) are 


9-7T 

-S L w' - d x K X = , 

Re 

_ Vl = -Hgw . 

R.iPe^ + 7 t 3 — Slu' — d z rci = —r^\7 2 w' , 


Re 
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(59) 

(60) 
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We see that the multiplier 7Ti plays a role similar to pressure, a standard result. Comparing 


the fourth equation with the constraint (13) also reveals that 7r 3 = RiPe0, which implies 


that we can eliminate both 9 and 7r 3 to get 


2RiPeV~V - S L u' - d z 7n = 


2tt- 


Re 


2 vV. 


(63) 


Dotting equations (59) to (61) with u' and integrating over the domain, we get (using 


incompressibility) the relationship: 


(, S L u'w') + RiPe^'V 2 ^) = ^-(p X7u'\ 2 ) = tt 2 V , 

Re 


(64) 


which reveals the interpretation of 7T2, and allows us to write 


T^lpn [ u ’} — (tt2 — 1)^ • 


(65) 


Since V is positive, this expression shows that energy stability corresponds to 7 t 2 < 1. All 


that is left to do is to solve equations (59)-(62) as well as the incompressibility condition for 













the eigenvalue 7T2 and determine for which values of Re and RiPe the condition 712 < 1 is 
satisfied. Unfortunately, this system of equations does not generally lend itself to Squire’s 
transformation, which means that we need to study the full 3D eigenproblem to solve for 
u', v', w', 9 , 7Ti and of course 7t 2 . 



FIG. 4. Stability boundaries for the LPN equations. The linear stability boundary is valid for 
infinite domains in both 2D and 3D. The energy stability boundary is shown in 2D (green curve) 
for a domain of arbitrary horizontal extension and in 3D (red curve) for the periodic domain of 
size L x x L y x L z = 10ir x 2ir x 2ir used for the low-Pe numerical simulations. The 3D energy 


stability limit falls between the lower and upper bounds (46) and (56). At large Reynolds number, 


the flow is linearly stable above a critical value (RiPe)^, and the energy stability limit follows the 
scaling (RiPe )e ~ Re. The symbols mark the simulations for which a turbulent solution was found 
numerically, using the LPN equations. 

Since Sl{z) = cos(z) for the Kolmogorov flow studied here, we use Floquet theory again 


to solve (59)-(62), together with -Ki = p ', 1 T 3 = RiPe# and the incompressibility condition. 


For simplicity, and for ease of comparison with the numerical simulations of the next section, 
we now restrict our attention to domains with vertical extent L z = 2n by setting the Floquet 
coefficient a = 0. Assuming an ansatz of the form 
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_ ilx-\-imy 


n=—N 


( 66 ) 











for each of the unknown functions yields the system: 



(69) 

(70) 

(71) 


( 68 ) 


( 67 ) 


w n + ( l 2 + m 2 + n 2 ) 6 n = 0 , 


lu n + mv n + nw n = 0 , 


for n = —N..N. This forms a generalized eigenvalue problem Az = 7 t 2 _Bz, where 
z = (u_ N , •••, u N ,v_ N , ....v N , w_ N , ..., w n ,P-n, @-n, ■■■, @ N ), which can be solved nu¬ 


merically (using LAPACK routines) for the eigenvalue 7r 2 . The latter depends on the hor¬ 


izontal wavenumbers l and m as well the original parameters Re and RiPe. For given Re 
and RiPe, energy stability is achieved if 7r 2 < 1 for all possible l and m. At fixed RiPe, 
the critical Reynolds number for energy stability is the largest value of Re for which this 
is true. Conversely, at fixed Reynolds number, the critical Richardson-Peclet number for 
energy stability (RiPe)# is the smallest value of RiPe for which 7r 2 < 1. The energy stability 
boundary therefore delimits the region of parameter space where the maximum value of 7 t 2 
over all possible l and m is smaller than unity. 

Figure [4ji compares the linear stability boundary to the energy stability boundary. The 
former is computed for an infinite domain and is valid both in 2D and in 3D, while the 
latter is computed for a 2D (y— independent) domain of infinite horizontal extent and for 
the 3D domain of size L x x L y x L z = 107r x 2 tt x 27r used in the low Peclet numerical 
simulations of the next section. Note how the 3D energy stability curve is lower than the 
2D one, which is expected since the family of all possible 2D perturbations is included 
in the family of all possible 3D perturbations. Systems whose parameters lie below the 
3D energy stability curve are always stable to perturbations of arbitrary amplitude, while 
systems whose parameters lie above the linear stability curve are unstable to infinitesimal 
perturbations. The 3D energy stability curve lies strictly below the linear instability curve, 
revealing a significant region of parameter space between them where a stratified shear flow 
is stable to small perturbations but could be destabilized by appropriate finite-amplitude 
perturbations. 

At large Reynolds number, (RiPe) e scales as Re with a proportionality constant of order 


unity, in agreement with the predictions from the upper and lower bounds. It is interesting 
to note that (RiPe)^; ~ Re is equivalent to (RiPr)# ~ 0(1), which is reminiscent of the 
nonlinear stability criterion originally proposed by Zahn®, albeit with the right-hand-side 
constant of order unity rather than of order 10~ 3 . His original argument, modified to have 
Re cr it = 1, could therefore provide a plausible physical explanation for the energy-stability 
scaling found. 

In summary, we have formally proved, using both simple analytical bounding arguments 
and exact numerical integration of the Euler-Lagrange equations, that a strong enough 
stratification makes the laminar shear flow stable to perturbations of arbitrary form and 
amplitude, within the constraint that the perturbations must still have a low Peclet number 


(see discussion in Section IVB2). 


V. NUMERICAL SIMULATIONS 

Our findings strongly suggest that stratified shear flows are subject to finite-amplitude 
instabilities, which raises the question of the relevance of linear stability analyses in de¬ 
termining when turbulent mixing is expected. In order to clearly assess the existence of 
finite-amplitude instabilities, we now turn to direct numerical simulations. 


A. The numerical model 


We solve the set of equations @-0i in a triply-periodic domain of size L x = ICR, L y = 47T 
and L z = 2i r, using the pseudo-spectral code originally developed by S. Stellmach to study 
double-diffusive convection 40 11 . This code has been modified to include the effect of the 
body force F. Table 1 shows a record of simulations run in this format. In all of these runs, 
Re = 10 4 , and Pe is either 0.1, 1 or 10. 


We then modified the code to solve instead the LPN momentum equation (16) together 


with the continuity equation, and have run a number of simulations with this new setup in 
a somewhat smaller domain ( L x = 107r, L y = 27T and L z = 2i r), for Re ranging from 10 2 to 
10 4 . The difference in the two domain sizes does not appear to have any influence on the 
numerical results in the low-Peclet-number regime, hence our decision to save on computer 
time in this second set of runs. The latter are summarized in Table 2. 






Pe 

Ri 

RiPe 

Transition to turbulence 

0.1 

1 

0.1 

Linear Instab. 

0.1 

10 

1 

Fin. amp. instab, starting from Ri = 1 run. 

0.1 

12 

1.2 

Fin. amp. instab, starting from Ri = 10 run. 

0.1 

15 

1.5 

No instab, found starting from Ri = 12 run. 

1 

0.001 

0.001 

Linear Instab. 

1 

0.01 

0.01 

Linear Instab. 

1 

0.1 

0.1 

Linear Instab. 

1 

0.3 

0.3 

Fin. amp. instab, starting from Ri = 0.1 run. 

1 

0.5 

0.5 

Fin. amp. instab, starting from Ri = 0.3 run. 

1 

0.7 

0.7 

Fin. amp. instab, starting from Ri = 0.5 run. 

1 

1 

1 

Fin. amp. instab, starting from Ri = 0.7 run. 

1 

1.2 

1.2 

Fin. amp. instab, starting from Ri = 1 run. 

1 

1.5 

1.5 

No instab, found starting from Ri = 1.2 run. 

10 

0.0001 0.001 

Linear Instab. 

10 

0.001 

0.01 

Linear Instab. 

10 

0.01 

0.1 

Linear Instab. 

10 

0.1 

1 

Fin. amp. instab, starting from Ri = 0.01 run. 

10 

0.12 

1.2 

Fin. amp. instab, starting from Ri = 0.1 run. 

10 

0.15 

1.5 

No instab, found starting from Ri = 0.12 run. 


TABLE I. Presentation of the various runs performed using the standard equations. All runs are at 
Re = 10 4 , in rectangular domains of size 107T x 47t x 27r. The resolution (in terms of equivalent mesh- 
points N x ,y, N z ) is the same in all directions, and for all runs, and is equal to 192 mesh points per 
interval of length 27r. Runs that go unstable starting from infinitesimal perturbations are marked 
“Linear Instab.”. Runs that do not go unstable starting from infinitesimal perturbations, but that 
have finite-amplitude instabilities are marked “Fin. amp. instab.”. These runs were started using 
the endpoint of another simulation at lower Ri, also noted. 



Re 

RiPe 

Transition to turbulence 


100 

0.1 

Linear Instab. 


100 

0.2 

Linear Instab. 


100 

0.22 

Linear Instab. 


100 

0.25 

No instab, found starting from RiPe = 

0.22 

1100 

0.09 

Linear Instab. 


1100 

0.21 

Linear Instab. 


1100 

0.24 

Linear Instab. 


1100 

0.27 

No instab, found starting from RiPe = 

0.24 

2500 

0.06 

Linear Instab. 


2500 

0.2 

Linear Instab. 



0.2 run. 
0.3 run. 
0.4 run. 
0.5 run. 
0.6 run. 
0.7 run. 
0.8 run 


2500 0.3 Fin. amp. instab, starting from RiPe = 

2500 0.4 Fin. amp. instab, starting from RiPe = 

2500 0.5 Fin. amp. instab, starting from RiPe = 

2500 0.6 Fin. amp. instab, starting from RiPe = 

2500 0.7 Fin. amp. instab, starting from RiPe = 

2500 0.8 Fin. amp. instab, starting from RiPe = 

2500 0.9 No instab, found starting from RiPe = 

Linear Instab. 

Linear Instab. 

tab. starting from RiPe = 0.1 run. 
tab. starting from RiPe = 0.3 run. 
tab. starting from RiPe = 0.5 run. 
tab. starting from RiPe = 0.6 run. 
tab. starting from RiPe = 0.8 run. 
stab, starting from RiPe = 1 run. 
und starting from RiPe =1.2 run. 


10000 

0.01 


10000 

0.1 


10000 

0.3 

Fin. amp. 

10000 

0.5 

Fin. amp. 

10000 

0.6 

Fin. amp. 

10000 

0.8 

Fin. amp. 

10000 

1 

Fin. amp. 

10000 

1.2 

Fin. amp 

10000 

1.5 

No instab. 


TABLE II. Presentation of the various runs performed using the LPN equations. All runs are in 
rectangular domains of size 107T x 27 t x 27 t. Those with Re = 10 4 have the same effective resolution 
as in Table 1. Those with Re = 100 have 96 meshpoints per interval of length 2ir and those with 
Re = 1100 and 2500 have 144 meshpoints per interval of length 27 t. 



B. Typical results 


We first take a look at typical simulations run using the LPN equations. Figure [5] shows 



FIG. 5. Snapshot of the streamwise (left) and vertical (right) velocity components for the run 
at Re = 10 4 and RiPe = 0.01 using the LPN equations, taken once it has equilibrated into a 
statistically-steady turbulent state. 

a system snapshot of a run at Re = 10 4 and RiPe = 0.01, once it has equilibrated into 
a statistically-steady turbulent state. The shear flow is visible in the left panel (which 
shows the velocity held in the x— direction), and the typical size and amplitude of the 
velocity perturbations are illustrated in the right panel (which shows the velocity held in the 
z—direction). For this particular value of RiPe, the shear is linearly unstable. We find that 
whenever this is the case, the system eventually settles into a statistically-steady turbulent 
state that is independent of the initial conditions. This is demonstrated in Figure [6^, which 
shows ( w 2 ) as a function of time for two simulations at Re = 10 4 and RiPe = 10 -4 : one 
that was started from small amplitude random initial conditions, and one that was started 
from the statistically-steady state reached by a previous run at Re = 10 4 and RiPe = 0.01. 
In both cases, ( w 2 ) settles into the same statistically-steady state after a short transient 
period. The same statement applies to all global diagnostics of the system dynamics. 

As shown in Figure [4^,, for Re = 10 4 the largest value of RiPe for which the laminar 
steady state solution ul(z) is linearly unstable is roughly equal to (RiPe)^ = 0.25. We 
have found that all low-Peclet-number simulations (i.e those run using the LPN equations, 
and those run with the standard equations at Pe < 1) which have RiPe < 0.25 do indeed 
transition to a turbulent state, and the results shown in Figure [5] and Figure^ are fairly 
representative of their behavior. A detailed quantitative analysis of the results of these runs 


will be presented elsewhere. 



FIG. 6. Evolution of ( w 2 ) as a function of time, for Re = 10 4 and RiPe = 10~ 4 (left) and 
RiPe = 0.5 (right) starting from small random perturbations (red solid line) and from finite- 
amplitude perturbations (green dashed line) by continuation of a previous run at other parameter 
values as noted in the legend. 

We have also found that this body-forced stratified shear flow is subject to hnite- 
amplitude instabilities for (RiPe)^ < RiPe < (RiPe) c , where the critical value (RiPe) c 
is discussed in the next section. This is shown in Figure [6}}, which presents ( w 2 ) as a func¬ 
tion of time for two simulations at Re = 10 4 and RiPe = 0.5: one that was started from weak 
amplitude random initial conditions, and one that was started from the statistically-steady 
state reached by a previous run at Re = 10 4 and RiPe = 0.1. We clearly see that the energy 
in the perturbations decays in the first case, but reaches a different statistically-steady state 
in the second case, a classical example of finite-amplitude instability. 

C. Finite-amplitude instability 

We now consider both the standard equations at Pe = 0.1, Pe = 1 and Pe = 10 and the 
LPN equations. In order to find turbulent solutions for RiPe > (RiPe )l more systematically, 
and determine the critical value (RiPe) c for the existence of finite-amplitude instabilities, 
we gradually increase Ri (keeping all other parameters fixed), using as a starting solution 
the result of a simulation run at lower Ri. We say that (RiPe) c is reached when we are no 
longer able to continue increasing Ri without losing the turbulent solution. Note that this 









only yields a rough estimate of (RiPe) c that depends largely on the size of the increments 
in Ri taken. It is possible that by using smaller increments one may be able to push further 
into the linearly stable region. Unfortunately, this is computationally very demanding and 
the increment size is in practice selected to satisfy our constraints on computation time. 

The results are summarized in Tables 1 and 2 and in Figure [4j and raise a number of 
interesting points. First note how (RiPe) c is the same for the LPN equations and for the 
standard equations at Pe = 0.1, Pe = 1 and even for Pe = 10 for Re = 10 4 . In all cases, 
we have (RiPe) c ~ 1.2. This validates the use of the LPN equations as a substitute for 
the standard equations for low-Peclet-number systems. One may in fact be surprised at the 
fact that the LPN equations even appear to be a good approximation of the large Pe runs 
(Pe = 10 here). However, this is due to the fact that the global Peclet number based on the 
amplitude of the hypothetical laminar shear flow ul is not a good predictor for the actual 
Peclet number of the turbulent solutions, Pe rms (see Section IIC). The latter is significantly 
smaller, and remains below one in all runs at Pe = 10. This result is consistent with the 
theory of LignieresOS, which merely requires Pe rms -C 1 for the the LPN equations to be 
valid. 

The value of (RiPe) c ~ 1.2 found at Re = 10 4 is somewhat larger than the linear stabil¬ 
ity threshold (RiPe)^, but is significantly smaller than the theoretical energy stability limit 
(RiPe)^ found in Section IV (see Figure [4]). Some level of discrepancy is expected, since 
lying within the energy stability limit is only a necessary (but not sufficient) condition for 
instability: even though everywhere within the energy-unstable domain there exist pertur¬ 
bations whose amplitude initially increase with time, this does not guarantee the onset of 
turbulence, as in most cases transient growth is followed by rapid decay. 

These results, however, show that at Re = 10 4 neither linear stability nor energy stability 
thresholds are good estimates for the actual threshold for transition to turbulence. Varying 
the Reynolds number from 100 to 10,000 we found that this is not always the case: for 
Re = 100 and Re = 1100, (RiPe) c and (RiPe)^ do appear to coincide and no finite-amplitude 
instabilities were found. The latter only appear for Re = 2500, and seem to exist at this 
Reynolds number for (RiPe)^ ~ 0.25 < RiPe < (RiPe) c ~ 0.8. 

The very limited finite-amplitude data available is not inconsistent with (RiPe) c r\_/ 0(1) 
for Re > 2500. We therefore see that, when using a non-dimensionalization based on the 
velocity and scale of the laminar flow, both the linear stability limit and the threshold to 




finite-amplitude instabilities are independent of the Reynolds number for large Re (at least, 
tentatively for the finite-amplitude threshold). The latter extends somewhat the stability 
threshold from the linear one (RiPe)^ ~ 0.25 to (RiPe) c ~ 0(1), but not by a large amount. 

VI. SUMMARY AND CONCLUSION 

We have analyzed in this work the stability of an idealized stratified, body-forced, low- 
Peclet-number shear flow using three different techniques: linear stability analysis, energy 
stability analysis, and direct numerical simulations. Our mathematical goal was three¬ 
fold: to test the validity of the LPN equations proposed by Lignieres®, to determine the 
respective thresholds for linear instability and energy stability, and to characterize the region 
of parameter space where finite-amplitude instabilities exist. 

Using dimensionless numbers based on the typical velocity of the laminar solution, our 
linear stability analysis confirmed that the LPN equations are indeed an excellent approx¬ 
imation to the standard equations of fluid dynamics provided Pe is smaller than 1. The 
domain of validity of these equations is in fact somewhat larger, and depends more on the 
Peclet number of the realized turbulent flow than the one of the hypothetical laminar so¬ 
lution. In the low Peclet number limit, thermal diffusion acts to destabilize the flow. We 
have found, as first shown by Lignieres' l4: and Lignieres, Califano, and Mangeney^, that 
the relevant bifurcation parameter is the Richardson number times the Peclet number, with 
stability for large Reynolds number achieved whenever RiPe > (RiPe)^ ~ 0.25. This shows 
that shear instabilities can exist at relatively large Richardson numbers in the small Peclet 
number limit. We have also shown using an extension of Squire’s transformation that in the 
same limit the first modes to be destabilized are 2D modes, a result which by contrast is 
not necessarily true for high-Peclet-number flows. 

We then performed an energy stability analysis of the LPN equations. We proved 
rigorously that any smooth low-Peclet-number shear flow becomes energy stable above a 
(Reynolds dependent) critical intensity of the background stratification. This has funda¬ 
mental implications: in this region of parameter space, the laminar flow is the only attractor 
of the dynamics, and therefore sustained turbulent mixing cannot take place. The criterion 
for energy stability is approximately RiPe > Re for large Reynolds number. Hence a lami¬ 
nar flow subject to strong stratification (with Ri > Pr -1 ) is energy stable, and the vertical 


diffusion of a scalar is due to molecular diffusivity only. 


These linear stability and energy stability results, however, may only be of academic 
interest. Indeed, using direct numerical simulations we have found that finite-amplitude in¬ 
stabilities in these low-Peclet-number stratified shear flows exist, for large enough Reynolds 
number, beyond the threshold for linear instabilities (RiPe)i ~ 0.25, but nevertheless disap¬ 
pear for RiPe significantly below the threshold for energy stability (RiPe )e ~ Re. Our very 
limited data is consistent with a finite-amplitude instability threshold (RiPe) c ~ 0(1) for 
large enough Re, using a non-dimensionalization based on the laminar velocity. These scal¬ 
ing laws are very tentative, in the sense that much remains to be done to measure (RiPe) c for 
larger Reynolds number, and to confirm the values found here. Indeed, as discussed above, 
it is possible that with more appropriately chosen initial conditions, one may be able to find 
turbulent solutions for even larger RiPe for a given Re. Furthermore, we note that while 
including the effects of rotation will not change the results of the energy stability analysis, 
it might allow for a wider range of dynamics and could help maintain turbulent solutions for 
larger RiPe. On the other hand, it is also not impossible that rotation could instead reduce 
the instability domain, or that some of the turbulent solutions found far into the region of 
linear stability are long chaotic transients that would eventually settle back to the laminar 
state upon longer numerical integration. Since the numerical constraints on the timestep 
and resolution increase dramatically for large Re and large RiPe flows, the accurate and 
definitive determination of (RiPe) c at large Reynolds number is a formidable task, one that 
should nevertheless be undertaken in the future. 

Finally, it is also worth recalling that all of these results only apply to the low Peclet 
number regime. While sufficiently-small-scale stellar shear layers fall into that category, 
large-scale shear layers, however, commonly have a high Peclet number (albeit still with a 
small Prandtl number). Both linear and energy stability analyses remain to be done in this 
case, and may reveal further surprises. 
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